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'. ABSTRACT 

^ ■ Classically, optical and near-infrared interferometry have relied on closure phase techniques to produce images. 
0^ [ Such techniques allow us to achieve modest dynamic ranges. 

I In order to test the feasibility of next generation optical interferometers in the context of the VLTI-spectro- 

■ imager (VSI), we have embarked on a study of image reconstruction and analysis. Our main aim was to test the 

, influence of the number of telescopes, observing nights and distribution of the visibility points on the quality of 

' the reconstructed images. Our results show that observations using six Auxiliary Telescopes (ATs) during one 

OO ! complete night yield the best results in general and is critical in most science cases; the number of telescopes is 

' the determining factor in the image reconstruction outcome. 

^ ' In terms of imaging capabilities, an optical, six telescope VLTI-type configuration and ~200 meter baseline 

k>( I will achieve 4 mas spatial resolution, which is comparable to ALMA and almost 50 times better than JWST will 

5^ ' achieve at 2.2 microns. Our results show that such an instrument will be capable of imaging, with unprecedented 

5^ I detail, a plethora of sources, ranging from complex stellar surfaces to microlensing events. 
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1. INTRODUCTION 



Image reconstruction is a key problem in observational astronomy, in particular in optical interferometry, where 
it is needed to overcome existing ambiguities in the interpretation of simple visibility measurements. 

In an array of TV telescopes, signals axe combined in ^ x N x (N-1) pairs or baselines to obtain ^ x N x 
(N-1) measurements called complex visibilities. These visibilities are related to the object brightness distribution 
via the van Cittert-Zernike theorem: 

V{u, v) = J J I{x, y) exp[—2 tt i {ux + vy)]dx dy 

where x and y arc angular displacements on the plane of the sky with the phase center as origin, I(x,y) is the 
brightness distribution of the target and u and v are the position vectors of the baselines projected on a plane 
perpendicular to the source direction, which together define the uv plane. In practical terms, the better the 
sampling of the uv plane in terms of baseline length, position angle, and number of measurements, the more 
faithful the reconstructed image will be relative to the true brightness distribution. 

In order to recover the brightness distribution, several image reconstruction algorithms for interferometry 
have been developed: phase referencing using AiPS, hybrid mapping using difmap (Masoni et al. 2005, Masoni 
2006), the building block bispectrum method (Hofmann & Weigelt 1993), mira (Thiebaut 2005), bsmem (Baron 
2007; private communication), among others. 

Closure phase imaging is the standard method for imaging in optical interferometry. It has been applied since 

the 1980's, initially to speckle interferometry (Weigelt ct al. 1998) and masked instrumentation and lately to 
optical interferometry with separate telescopes. This technique relies on two quantities related to the complex 
visibility - the squared visibility and the triple product - to perform image reconstruction. 

The squared visibility or power-spectrum is defined as: 

V{ui,Vi)V*{ui,Vi) = \V{ui,Vi)\'^ 

and the triple product or bi-spectrum: 

T{ui,Vi,Uj,Vj) = V{ut,Vt)V{uj,Vj)V*{ui +Uj,Vi + Vj) 

where Ui+Uj=Uk, Vi+Vj=Vk, Ui, Vj, % ,Vj, are vectors defining the position of three baselines and the asterisk 
denotes the complex conjugate. The triple product is the product of the complex visibilities on the baselines 
forming a closed loop joining any three telescopes (Fig. 1). The modulus of the triple product is the triple product 
amplitude and the argument is the closure phase. The closure phase is independent of atmospheric-induced phase 
variations. 




Figure 1. 



2. ARRAY CONFIGURATION 



AMBER was the first VLTI instrument to use beam recombination for three telescopes. Next generation instru- 
ments hke VSI (Malbet et aL 2006) and Matisse (Lopez et ah 2006) wiU use a larger number of baselines. In 
order to test the importance of the number, type and configuration of the telescopes and number of observing 
nights on the image reconstruction, we have adopted three telescope setups using the VLTI as a template. The 
setups were chosen to optimize the uv plane coverage, while taking into consideration such limitations as the 
delay lines. 

The following configurations were adopted: 




Figure 2. 4 UT x 1 niglit configuration. 
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Figure 3. 4 AT x 3 night configuration (left). 6 AT x 1 night configuration (right). 



• 4 UTs x 1 night (U1-U2-U3-U4) - faint source configuration (Fig. 2); 

• 4 ATs X 3 nights (A0-B1-D2-G1, Dl-EO-HO-JB, AO-GO-Il-MO) - the minimum imaging uv coverage con- 
sistent with the current VLTI capabilities (Fig. 3; left); 

• 6 ATs X 1 night (A0-B1-D2-G1-J2-M0) - a 6 telescope extended configuration with one night observation; 
has less uv points than the 4 AT x 3 nights configuration (Fig. 3; right). 



3. NOISE MODEL 



The noise model estimates the uncertainties on squared visibilities and closure phases assuming an instrument 
following a multi-axial recombination scheme with a fringe tracker (FT; Jocou et al.2007). 



The quantity we procure is the total number of detected photoevents per basehne in 5 minutes in the 
photometric channel: 



)V _ ^tot»tx(l-/)x5xJVofc, 

where / is the fraction of the beam that gores into the intcrfcrometric channel (90%), Nobs is the number of 
observations per minute, Nbase is the number of baselines and Ntotai is the total number of detected photevents 
per integration point in all channels: 

Ntotai = FqX iQ-O-iniag ^ ^.^^ ^ jy^^^ X TT X i?^ X AA X traus X Strehl 

Here _Fo is the photon flux of a zcroth magnitude star (Jocou ct al. 2007), mag is the object magnitude in the 
observing band, Unt is the integration time (0.05 s), Nt^i is the number of telescopes (6 or 4), i? is the radius 
of the telescopes (4.1 for UTs and 0.9 for ATs) assumed for simplification to have no central hole, AA is the 
spectral bandwidth (chosen), trans is the total instrument transmission including quantum efficiency (Jocou et 
al. 2007), and Strehl is the Strehl ratio and depends on wavelength (Jocou et al. 2007). 

The correlated flux per baseline is given by: 

^ cor — 2 

where V is the instrumental visibility loss (80%) and FT loss (90%) corrected visibility and Ni = Ntotai x / x 
5 X Nobs /Nbase, is the detected number of photoevents per baseline in 5 minutes in the interferometric channel. 

The signal-to- noise ratio of the squared correlated flux (SNR) is a complex function of Fcor, Ni, V , Ntei, 
Npix the number of pixels needed to code the interferometric channel (600) and a the detector readout noise (15 
e-). 

Therefore, the uncertainty on the squared visibility is: 

errorv2 = V'^ x ^SNR-^ -h 2 x (iVp + N^i^p a^/N^) 

where Npi^p is the number of pixels needed to code the photometric channel (4). 

The photon noise on the closure phase {Phot) is a function of Nteu Ni, V123, Vi, V2, V3, V12, ¥2:^, Vis, where 
the V terms correspond to visibilities in a combination of one, two or three telescopes. The detector noise on 
the closure phase (Det) depends on Ntei, Ni, Npi^, a and the V terms. 

Therefore, the total noise on the closure phase is estimated as: 

4. OIFITS FILE GENERATION 

The advent of optical interferometry has lead to the need for a common data format based on the Flexible 
Image Transport System (FITS) in order to facilitate data reduction and the combination of various types of 
interferometer data. The result was the creation of the 01 Exchange Format or OIFITS format (Pauls et al. 
2004). Interferometer projects supporting this standard include COAST, NPOL, IOTA, VLTI, PTI and the 
Keck Interferometer. 

OIFITS format is designed so that different categories of information arc stored in distinct tables within a 
file and can be cross-referenced one to another. Each OIFITS file table stores specific parameters that include 
interferometric observables - complex visibility (01- VIS), squared visibilities (0I-VIS2), and the triple product 
(bispectrum/closurc phase; 0I-T3) - and ancillary data like Universal Time measurements, spectral bandpass and 
wavelengths (01- WAVELENGTH), target identification (OI-TARGET) and telescope coordinates (01- ARRAY). 



Synthetic images were provided by the science case groups (Garcia et al. 2007) and imported into ASPRO, an 
image simulation tool originally created for IRAM. Assuming the objects were located at a —60 degree declination 
angle, typical object magnitudes were used to generate noisy, squared visibilities and closure phases. We assumed 
an one hour per calibrated measurement. The actual on-source integration time is, however, 10-15 minutes per 
hour due to overheads. The total integration time assumes an entire transit (9 hours). 

5. PHASE CLOSURE THEORY 



1- 

V 

Figure 4. 

Given an incoherent source, the observed brightness distribution in a direction s' is given by (Fig. 4): 

jobs ^ jobs (^x,y) ^ PSF{x, y) * /*™^ {x,y)+N{x,y) 

where PSF(x,y) is the point spread function, f"^^^ (x,y) is the true object brightness distribution, N(x,y) is the 
noise and the asterisk denotes convolution. 

In practice, interferometry does not make measurements in the image plane but in Fourier space. The relevant 
quantity is called the complex visibility and is measured at each uv point, the position vector of the baseline on 
a plane perpendicular to the source direction: 

V^'^iu^v) = y*™'=(u,w) X S{u,v) + N'{u,v) 

where S(u,v), the sampling function, is the Fourier transform of the PSF(x,y), \^^"'^(x,y), the true visibility, is 
the Fourier transform of the true brightness distribution f™'' (x,y) and N'(u,v) is the noise in the Fourier space. 

The van Cittert-Zernike theorem states that the true brightness distribution can be obtained by the inverse 
Fourier transform and deconvolution of the observables: 

P'^''%x,y)*PSF{x,y) = / / y*''"'=(u.w) x S{u,v) exp[2ni{ux + vy)] dudv 

The role of image reconstruction is to obtain the best approximation, P'P'^°^ (x,y) ~ l*™^(x,y), to the true 
brightness distribution. For phase closure image reconstruction we want to find a solution which has the maximum 
probability given the data. The solution P'P'^°^ (x,y) involves minimizing the value of a penalty term subject to 
some condition: 

P{z) = Pl{z) + Pp„or{z) = Pl{z) + tiR{z) 

where z represents the intensity of the pixels, Pl is the likelihood penalty, Pprior is the prior penalty, R(z) is the 
regularization term and /i is a multiplier (hyperparameter) tuned so that at the solution, the likelihood terms are 
equal to their expected values. Pl enforces agreement with the data, while Pprior provides information where 
the data fail to do so, in particular in regions where the noise dominates or where the data are missing. Pprior 
is responsible for the so-called regularization of the inverse problem. 

The solution is therefore found by minimizing the likelihood penalty term for squared visibilities and triple 
products under a prior constraint: 



where is the likelihood term with respect to the squared visibility data and to the triple product. 
Typically, the data penalties are defined assuming the measurements follow Gaussian statistics: 

Xv^ (^) = — 2~ (Klata ~" ^model) 



where refers to the squared visibilities (data and model), (/)^"*" to the measured closure phase, tpi, ip2, <^3 to 
three individual phases involved in the phase closure and a to the respective standard deviations. 

Regularization can enforce agreement with some preferred and/or exact properties of the solution. Because 
we want to favor compact sources, smoothness is a very common regularization constraint: 

R{z) = Y.w\z\'' 

where the regularization weights w are chosen in order to ensure spectral smoothing, that is, it enforces smooth- 
ness of the Fourier spectrum of the image and enforces compactness of the brightness distribution of the field of 
view. 

We have implemented the Multi-Aperture Reconstruction Algorithm (mira; version 0.7; April 2008) as 
developed by Eric Thiebaut with support from the Jean-Marie Mariotti Center (JMMC). MIRA works in the 
YORICK and C platform and is optimized to handle optical interferometric data with sparse uv coverage (Thiebaut 
2005). Typically, MIRA implements several types of regularization of which smoothness and maximum entropy 
are some examples. 

6. IMAGE ANALYSIS 

Because mira typically superresolves in the image reconstruction, we have compared our reconstructed images 
with the synthetic images provided by the science groups. Relative astrometry information was obtained for the 
images using DS9. Photometry of the image components was performed using IRAF procedure phot. SNRy2 and 
SNRt correspond to the range of signal-to-noise ratios in the squared visbilities and closure phases measured 
in the synthetic images. 

7. CLOSURE PHASE IMAGE RECONSTRUCTION 

Table 1. AGN. Diameter units are in pixels. 





Image 


MIRA 4 UT 


flux nucleus 


22.3% 




flux inner diameter 


55.0% 




flux outer diameter 


17.6% 


91.6% 


flux torus 


5.1% 


8.4% 


nuclear diameter 


7 




inner diameter 


30 




outer diameter 


75 


70 


torus diameter 


260 




SNRv2 




553-704 


SNRt 




0.3-0.4 



milliarcsec 



milliarcsec 



Figure 5. A simulated AGN with torus, 0.1 mas/pixel sampling; mira reconstruction 4 UT x 1 night configuration, 
smoothness regularization, 0.1 mas/pixel sampling. 
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milliarcsec 



Figure 6. A low mass evolved star system with an outflow, 0.1 mas/pixel sampling; mira reconstruction 4 AT x 3 nights 
configuration, smoothness regularization, 0.1 mas/pixel sampling; mira reconstruction 6 AT x 1 night configuration, 
smoothness regularization, 0.1 mas/pixel sampling. 



Table 2. Evolved Star. Diameter units are in pixels. 





Image 


MIRA 4 AT 3 


MIRA 6 AT 


fiux star 


25.1% 


22.3% 


22.5% 


fiux wind 


74.9% 


77.7% 


77.5% 


ratio star/ wind 


0.3 


0.3 


0.3 


star diameter 


10 


15 


12 


inner wind diameter 


60 X 40 


50 X 30 


60 X 40 


outer wind diameter 


100 X 80 


90 X 80 


100 X 80 


SNR^2 




63-2336 


38-1706 


SNRt 




0.01-0.2 


1.0-44 



E 




milliarcsec milliarcsec 



Figure 7. A simulated stellar surface of a M giant, 0.1 mas/pixel sampling; mira reconstruction 4 AT x 3 nights config- 
uration, smoothness regularization, 0.1 mas/pixel sampling, SNR„2 = 1. 5-541, SNRt=0. 01-05; mira reconstruction 6 AT 
X 1 night configuration, smoothness regularization, 0.1 mas/pixel sampling, SNRy2 =0.1-600, SNRt=0. 1-115. 
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Figure 8. A simulated stellar surface of a supergiant, 0.1 mas/pixel sampling; mira reconstruction 4 AT x 3 nights con- 
figuration, smoothness regularization, 0.1 mas/pixel sampling, SNR=y2 =4-4900, SNRt= 0.004-0.09; mira reconstruction 
6 AT X 1 night configuration, smoothness regularization, 0.1 mas/pixel sampling, SNRy2 =4-2999, SNRt=0. 2-116. 




Figure 9. A simulated microlensing event, 0.1 mas/pixel sampling; mira reconstruction 4 UT x 1 night configuration, 
smoothness regularization, 0.1 mas/pixel sampling. 

Table 3. Microlensing. Distance is in pixels. 





Image 


MIRA 4 UT 


A flux 


19.6% 


19.6% 


B fiux 


11.9% 


11.9% 


C flux 


29.8% 


29.8% 


D flux 


38.7% 


38.7% 


ratio D/A 


2.0 


2.0 


ratio D/B 


3.3 


3.3 


ratio D/C 


1.3 


1.3 


distance AB 


40 


40 


distance BC 


20 


20 


distance CD 


56 


56 


SNR^,2 




228-268 


SNRt 




0.1-0.4 



Table 4. Pinwheel 





Angle 


Image 


SNRv-2 


SNRt 


MIRA 6 AT 


inner spiral 
inner spiral 


deg 
60 deg 


100 x 80 
100 x 60 


49-861 
36-861 


2-19 
1-43 


100 X 80 
100 X 60 



8. DISCUSSION 

8.1 The Sources 

In general, the image reconstruction faired very well. This is denoted by the excellent astrometry and photometry 
of the individual image elements. Some sources merit a separate discussion. 

We are able to detect the unresolved AGN out to the sublimation radius. The total flux contained in 
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Figure 10. A simulated image of the pinwheel nebula at degree inclination, 1.0 mas/pixel sampling; mira reconstruction 6 
AT X 1 night configuration, smoothness regularization. A simulated iamge of the pinwheel nebula at 60 degree inclination, 
1.0 mas/pixel sampling; mira reconstruction 6 AT x 1 night configuration, smoothness regularization, 0.1 mas/prxel 
sampling. 
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Figure 11. A simulated stellar surface of the structure of inner disks surrounding YSOs, 0.1 mas/pixel sampling; mira 
reconstruction 4 UT x 1 night configuration, smoothness regularization, 0.1 mas/pixel sampling; mira reconstruction 4 
AT X 3 nights configuration, smoothness regularization, 0.1 mas/pixel sampling; mira reconstruction 6 AT x 1 night 
configuration, smoothness regularization, 0.1 mas/pixel sampling. 



Table 5. YSO. Radii units are in pixels. 





Image 


MIRA 4 UT 


MIRA 4 AT 3 


MIRA 6 AT 


flux star 


18.1% 


18.7% 


18.6% 


18.6% 


flux disk 


81.9% 


81.3% 


81.4% 


81.4% 


ratio 


0.2 


0.2 


0.2 


0.2 


inner diameter 


65 X 45 


65 X 45 


65 X 45 


65 X 45 


outer diameter 


85 X 60 


80 X 56 


82 X 60 


85 X 60 


SNRv2 




2055-5215 


89-1385 


12-1363 


SNRt 




0.04-0.2 


0.04-0.5 


0.2-12 



this region is 91.6% compared to the 91.4% measured in the synthetic image. The result suggests that such 
observations will provide important constraints on the size of the sublimation region and torus in the nearest 
AGN. 

For the stellar surfaces, our results show unprecedented detail. Most of the flux is recovered in a compact 
region and careful inspection of the reconstructed images shows correspondence with individual regions of different 
brightness in the synthetic image. The edges of the stellar surface are also well reproduced. 

8.2 Number of Telescopes and uv Coverage 

Our results show that the use of the 6 AT x 1 night configuration generally yields better reconstructed images 
than the 4 AT x 3 nights conflguration, even though the number of uv points is similar. This is particularly noted 
in images of stellar surfaces and in stellar disks/winds. This is in consonance with previous studies (Tuthill & 
Monnier 2000; Baron 2007). The image reconstruction result is satisfactory with a low number of array elements 
when the source is compact (Tuthill & Monnier 2000). However, when the sources are more complex, a larger 
Fourier plane coverage is needed and therefore more array elements are needed. 

Indeed it is the number of telescopes and not the number of observing nights (and therefore total number 
of uv points) that determine the image reconstruction outcome. In a system of A'^ telescopes, the number of uv 



points is X {N — 1) per night per integration point and the number of independent closure phases is | x 

(N-1) X (N-2) per night per integration point. With N=4 and 3 nights, for each integration point there arc a 
total of 18 uv points (therefore squared visibilities) and a total of 9 closure phases. For N=6 and 1 night, the 
number of squared visibilities and closure phases per integration point is 15 and 10, respectively. 

On the other hand, the image reconstruction (or fitting process) can only make use of the instantaneous phases 
to create a model brightness. Therefore, for a N=4 system, the image reconstruction uses the 3 instantaneous 
closure phases (3 per night), while with 6 telescopes there are now 10 independent closure phases. 

8.3 Dynamic Range 

Image reconstruction, whatever the method, is crucial when a high dynamic range (>1000) is required. Our 
results show that the closer phase technique allows to achieve modest dynamic ranges. Higher dynamic ranges 
can be obtained by using the VLTI coupled with a phase referencing instrument, like PRIMA (Dcplankc ct al. 
2000) which will allow to reach K band magnitudes of 11 mag for the ATs and 14 mag for the UTs. The presence 
of nearby phase referencing sources will allow larger integration times on source and therefore higher dynamic 
range. 

9. CONCLUSIONS 

Next generation optical interferometric instrumention like the VSI will be designed to offer superb imaging 
capabilities in order to bypass ambiguities in the interpretation of visibility data. It will be particularly novel in 
that it will allow spatial resolutions of milliarcseconds, comparable to ALMA and even better the future JWST. 

In order to test future VLTI performance in terms of imaging, we have run several tests on key scientific cases. 
Results show that with a simple 6 Auxiliary Telescope (AT) configuration and 1 night observation, images of high 
image fidelity can be obtained. Using the VLTI alone will allow to achieve modest dynamic ranges, although in 
conjunction with a phase referencing instrument, K band magnitudes of 11 and 14 mag can be reached with the 
ATs and UTs, respectively. Such diverse objects as AGN, stellar surfaces and disks are shown possible to image 
with future VLTI instumentation. 
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